# Replication package for 
# "The Economic Leverage of International Organizations in Interstate Disputes"
# Johannes Karreth
# June 30, 2017
# jkarreth@ursinus.edu

# This file: 23_crises_analyses.R
# Purpose: Analyses of major clashes or war during crises

# Note: must have read in data before, using 21_crises_readdata.R

# File structure:
# 1. Regression estimates, including some posterior predictive checks
# 2. Bayesian Model Averaging
# 3. IGO characteristics

# setwd("...")

library("rstan")
library("rstanarm")
library("texreg")
library("ggplot2")
library("xtable")
library("gridExtra")
library("loo")
library("lme4")
library("dplyr")
library("reshape2")
library("BMA")

# Source in functions
source("Functions/MCMClogit.fd.mat.R")
source("Functions/theme_jk.R")
source("Functions/plotBMA.R")
source("Functions/MCMC_roc_prc.r")
source("Functions/geom_flat_violin.R")

#############################################
# 1. Regression estimates, including some posterior predictive checks
#############################################

# Unit of analysis: crisis

# Model 1 (for Figure 2)

m1_eq <- viol.majclash ~ igo_lev3_count_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally
m1_mf <- model.frame(m1_eq, data = dat)

m1 <- glm(m1_eq, data = m1_mf, family = binomial(link = "logit"))

m1_stanglm <- stan_glm(formula = m1_eq,
                       data = (m1_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m1_stanglm, file = "Output_MCMC/crises_m1_stanglm.RDS")

# Table

m1_tab <- m1_stanglm$stan_summary[, c(1, 3)]
m1_tab <- data.frame(m1_tab[c(2:8, 1, 10), ])
m1_tab$varname <- c("IGOs with high leverage", 
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m1_tab) <- c("Mean", "Standard deviation", "Variable")
m1_tab <- m1_tab[, c(3, 1, 2)]

m1_xtab <- xtable(m1_tab, digits = 2)
print(m1_xtab, file = "Output_Tables-and-Figures/crises_m1_table.tex", include.rownames = FALSE)

# First differences (Figure 2)

m1_mm <- model.matrix(m1_eq, data = dat)

temp_sims <- as.matrix(m1_stanglm)
m1.fd <- MCMClogit.fd.mat(model_matrix = m1_mm, mcmc_out = temp_sims, 
                          credint = c(0.025, 0.975), 
                          percentiles = c(0.1, 0.9),
                          full_sims = TRUE)

m1.changes <- matrix(NA, ncol = 2, nrow = ncol(m1_mm))
for (i in 2:ncol(m1_mm)){
  m1.changes[i, ] <- ifelse(length(unique(m1_mm[, i])) == 2 & range(m1_mm[, i]) == c(0, 1), c(0, 1), 
                            quantile(m1_mm[, i], probs = c(0.1, 0.9)))
}

# > m1.changes
# [,1]  [,2]
# [1,]    NA    NA
# [2,] 0.000 5.000
# [3,] 0.000 1.000
# [4,] 0.000 1.000
# [5,] 0.000 1.000
# [6,] 0.000 1.000
# [7,] 0.143 4.055
# [8,] 0.000 1.000

fd.dat <- data.frame(m1.fd)
fd.dat$variable_label <- c("IGOs with high leverage\n(from 0 to 5)", 
                           "Existential threat\n(no to yes)", 
                           "Territorial dispute\n(no to yes)", 
                           "Joint democracy\n(no to yes)", 
                           "Strategic rivalry\n(no to yes)",
                           "UNGA ideal point difference\n(0.1 to 4)",
                           "Allies \n(no to yes)")

fd_long <- melt(m1.fd)
fd_long$Varpos <- ifelse(fd_long$Var2 == "igo_lev3_count_use", 1,
                         ifelse(fd_long$Var2 == "highgrav", 2, 
                                ifelse(fd_long$Var2 == "territory", 3,
                                       ifelse(fd_long$Var2 == "jointpol7", 4,
                                              ifelse(fd_long$Var2 == "rivalry_th", 5,
                                                     ifelse(fd_long$Var2 == "absidealdiff", 6, 
                                                            ifelse(fd_long$Var2 == "atopally", 7,
                                                                   NA)))))))

fd_long$variable_label <- ifelse(fd_long$Var2 == "igo_lev3_count_use", "IGOs with high leverage\n(from 0 to 5)",
                                 ifelse(fd_long$Var2 == "highgrav", "Existential threat\n(no to yes)", 
                                        ifelse(fd_long$Var2 == "territory", "Territorial dispute\n(no to yes)",
                                               ifelse(fd_long$Var2 == "jointpol7", "Joint democracy\n(no to yes)",
                                                      ifelse(fd_long$Var2 == "rivalry_th", "Strategic rivalry\n(no to yes)",
                                                             ifelse(fd_long$Var2 == "absidealdiff", "UNGA ideal point difference\n(0.1 to 4)", 
                                                                    ifelse(fd_long$Var2 == "atopally", "Allies \n(no to yes)",
                                                                           NA)))))))

# ROPE: use SE of the ratio of 1s

r_n <- length(m1_mf$viol.majclash)
r_p <- sum(m1_mf$viol.majclash) / r_n
r_se <- sqrt((r_p * (1 - r_p)) / r_n)
ROPE <- r_se

fd.p <- ggplot(data = fd_long, aes(x = reorder(variable_label, -Varpos), y = value)) + 
  scale_y_continuous(labels = function(x) x*100) +
  # scale_x_discrete(labels = rev(fd.dat$variable_label)) +
  # geom_rect(aes(ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf), col = "gray") +
  # geom_segment(y = 0, yend = 0, x = 0, xend = 10.5, color = "black", linetype = "solid", size = 0.33) + 
  geom_flat_violin(color = NA, fill = "black") +
  coord_flip() + 
  theme_classic() + 
  ylab("Percentage point change in the probability of major clashes\nduring a crisis as each predictor changes as indicated") + 
  xlab("") + 
  theme(axis.ticks = element_blank(), axis.text = element_text(color = "black"), axis.line = element_blank()) + 
  annotate(geom = "text", y = -0.05, x = 0.15, label = "Major clashes less likely", hjust = 1, size = 3) +
  annotate(geom = "text", y = 0.05, x = 0.15, label = "Major clashes more likely", hjust = 0, size = 3) +
  annotate(geom = "segment", y = 0, yend = -0.7, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "segment", y = 0, yend = 0.7, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "rect", ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf, fill = "gray", color = NA, alpha = 0.5) +
  geom_segment(y = 0, yend = 0, x = 0, xend = 10.5, color = "black", linetype = "solid", size = 0.33) + 
  annotate(geom = "text", y = -0.15, x = 5.25, hjust = 1, vjust = 0.5, label = "Region of practical\nequivalence (ROPE)\n[-2, 2]", color = "darkgray", size = 3) +
  annotate(geom = "segment", y = -0.15, yend = -1 * ROPE, x = 5.25, xend = 5.3, size = 0.25, color = "darkgray") +
  annotate(geom = "text", y = 0.1, x = 7.25, hjust = 0, vjust = 0.5, label = "% of posterior distribution\nto the left of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.1, yend = -0.2, x = 7.25, xend = 7.15, size = 0.25) +
  annotate(geom = "text", y = 0.35, x = 4.35, hjust = 0.5, vjust = 0, label = "% of posterior distribution\nto the right of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.3, yend = 0.2, x = 4.35, xend = 4.05, size = 0.25) +
  annotate(geom = "text", y = -0.45, x = 6.5, hjust = 0.5, vjust = 0, label = "Mean estimated change", size = 3) +
  annotate(geom = "segment", y = -0.5, yend = -0.3, x = 6.65, xend = 6.9, size = 0.25)

# % < 0

m1.fd_out0 <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0) / length(x), sum(x > 0) / length(x)))
m1.fd_outROPE <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0 - ROPE) / length(x), sum(x > 0 + ROPE) / length(x)))

fd_changelab <- ifelse(apply(m1.fd, 2, median) < 0, "Decrease by", "Increase by")
fd_changelab_short <- ifelse(apply(m1.fd, 2, median) < 0, "-", "+")
fd_changetext <- paste(fd_changelab, abs(round(apply(m1.fd, 2, median) * 100)), "points", sep = " ")
fd_changetext_short <- paste(fd_changelab_short, abs(round(apply(m1.fd, 2, median) * 100)), sep = "")

m1.fd_annotate <- data.frame(x = apply(m1.fd, 2, median), 
                             y = rev(unique(fd_long$Varpos)), 
                             out0 = paste(round(m1.fd_out0 * 100, digits = 1), "%", sep = ""), 
                             outROPE = paste(round(m1.fd_outROPE * 100, digits = 1), "%", sep = ""),
                             change = fd_changetext_short)
fd.p2 <- fd.p + geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = outROPE), color = "white", nudge_x = 0.075, size = 3) + 
  geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = change), color = "black", nudge_x = -0.1, size = 4)

ggsave(fd.p2, file = "Output_Tables-and-Figures/crises_m1_fd.pdf", width = 6, height = 8)

# Remove transparency
fd.p3 <- fd.p + annotate(geom = "rect", ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf, fill = "gray90", color = NA) + geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = outROPE), color = "white", nudge_x = 0.075, size = 3)  + 
  geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = change), color = "black", nudge_x = -0.1, size = 4)
ggsave(fd.p3, file = "Output_Tables-and-Figures/crises_m1_fd.eps", width = 6, height = 8)

# LOO

m1_noigo <- stan_glm(formula = viol.majclash ~ highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally,
                     data = (m1_mf),
                     family = binomial(link = "logit"),
                     prior = normal(0, 10), 
                     prior_intercept = normal(0, 10),
                     seed = 20170207,
                     cores = 4)

saveRDS(m1_noigo, file = "Output_MCMC/crises_m1_noigo_stanglm.RDS")

library("loo")
loo_noigo <- loo(m1_noigo, cores = 2)
loo_m1 <- loo(m1_stanglm, cores = 2)
compare(loo_noigo, loo_m1)
# elpd_diff se 
# 2.9       3.0

# Influence of individual observations
p <- ggplot(data = data.frame(x = 1:nrow(m1_mf), y = loo_m1$pareto_k),
            aes(x = x, y = y)) + geom_point(shape = 1) + xlab("Observation") + ylab("Shape parameter k") + geom_hline(yintercept = 0.5, linetype = "dashed") + theme_jk()
ggsave(p, file = "Output_Tables-and-Figures/crises_loo_k.pdf", width = 5, height = 4)

# Other model fit diagnostics via Stan

# Data 
roc_mf <-  model.frame(viol.majclash ~ igo_lev3_count_use + igo_ingr23_use + igo_mp3_use + igo_pb_use + igo_allothers + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally, data = dat)
N <- nrow(roc_mf)
majclash <- roc_mf$viol.majclash
igo_lev3_count_use <- roc_mf$igo_lev3_count_use
igo_ingr23_use <- roc_mf$igo_ingr23_use
igo_mp3_use <- roc_mf$igo_mp3_use
igo_pb_use <- roc_mf$igo_pb_use
igo_allothers <- roc_mf$igo_allothers
highgrav <- roc_mf$highgrav
territory <- roc_mf$territory
jointpol7 <- roc_mf$jointpol7
rivalry_th <- roc_mf$rivalry_th
absidealdiff <- roc_mf$absidealdiff
atopally <- roc_mf$atopally

m1_roc_data_list <- c("N", "majclash", "igo_lev3_count_use", "igo_ingr23_use", "igo_mp3_use", "igo_pb_use", "igo_allothers", "highgrav", "territory", "jointpol7", "rivalry_th", "absidealdiff", "atopally")

# Fit m1 directly in Stan
m1_roc_rstan <- stan(file = "crises_m1_roc.stan", 
                     data = m1_roc_data_list,
                     iter = 1000, chains = 4,
                     seed = 20170207,
                     cores = 4)
saveRDS(m1_roc_rstan, file = "Output_MCMC/crises_m1_rstan.RDS")

# Fit igo_ingr23_use directly in Stan
m_ingr23_roc_rstan <- stan(file = "crises_m-ingr23_roc.stan", 
                           data = m1_roc_data_list,
                           iter = 1000, chains = 4,
                           seed = 20170207,
                           cores = 4)
saveRDS(m_ingr23_roc_rstan, file = "Output_MCMC/crises_m1_ingr23_rstan.RDS")

# Fit igo_mp3_use directly in Stan
m_mp3_roc_rstan <- stan(file = "crises_m-mp3_roc.stan", 
                        data = m1_roc_data_list,
                        iter = 1000, chains = 4,
                        seed = 20170207,
                        cores = 4)
saveRDS(m_mp3_roc_rstan, file = "Output_MCMC/crises_m1_mp3_rstan.RDS")

# Fit igo_pb_use directly in Stan
m_pb_roc_rstan <- stan(file = "crises_m-pb_roc.stan", 
                       data = m1_roc_data_list,
                       iter = 1000, chains = 4,
                       seed = 20170207,
                       cores = 4)
saveRDS(m_pb_roc_rstan, file = "Output_MCMC/crises_m1_pb_rstan.RDS")

# Fit all other IGOs directly in Stan
m_allothers_roc_rstan <- stan(file = "crises_m-allothers_roc.stan", 
                              data = m1_roc_data_list,
                              iter = 1000, chains = 4,
                              seed = 20170207,
                              cores = 4)
saveRDS(m_allothers_roc_rstan, file = "Output_MCMC/crises_m1_allothers_rstan.RDS")

# Fit empty model (no IGOs) directly in Stan
m_noigos_roc_rstan <- stan(file = "crises_m-noigos_roc.stan", 
                           data = m1_roc_data_list,
                           iter = 1000, chains = 4,
                           seed = 20170207,
                           cores = 4)
saveRDS(m_noigos_roc_rstan, file = "Output_MCMC/crises_m1_noigos_rstan.RDS")

# ROC and precision-recall curves

m1_rp <- MCMC_roc_prc(stan_object = m1_roc_rstan, 
                      model_frame = m1_mf, 
                      model_name = "IGOs with high leverage")

m_ingr23_rp <- MCMC_roc_prc(stan_object = m_ingr23_roc_rstan, 
                            model_frame = m1_mf, 
                            model_name = "Structured IGOs")

m_mp3_rp <- MCMC_roc_prc(stan_object = m_mp3_roc_rstan, 
                         model_frame = m1_mf, 
                         model_name = "Highly structured Security IGOs")

m_pb_rp <- MCMC_roc_prc(stan_object = m_pb_roc_rstan, 
                        model_frame = m1_mf, 
                        model_name = "Peace-brokering IGOs")

m_allothers_rp <- MCMC_roc_prc(stan_object = m_allothers_roc_rstan, 
                               model_frame = m1_mf, 
                               model_name = "All other IGOs")

m_noigos_rp <- MCMC_roc_prc(stan_object = m_noigos_roc_rstan, 
                            model_frame = m1_mf, 
                            model_name = "No IGOs")

PRC_dat <- rbind(m1_rp$prc_dat, m_ingr23_rp$prc_dat, m_mp3_rp$prc_dat, m_pb_rp$prc_dat, m_allothers_rp$prc_dat, m_noigos_rp$prc_dat)
PRC_labeldat <- rbind(m1_rp$area_under_prc, m_ingr23_rp$area_under_prc, m_mp3_rp$area_under_prc, m_pb_rp$area_under_prc, m_allothers_rp$area_under_prc, m_noigos_rp$area_under_prc)

PRC.p <- ggplot(data = PRC_dat, aes(x = x, y = y)) + geom_line() + ylim(c(0, 1)) + facet_wrap(~ Model, ncol = 1) + geom_text(data = PRC_labeldat, aes(label = paste("Area under\ncurve:",round(auc, digits = 2))), x = 0, y = 0.2, hjust = 0) + labs(title = "Precision-recall curves", hjust = 0.5) + xlab("Recall") + ylab("Precision") + theme_jk()

ROC_dat <- rbind(m1_rp$roc_dat, m_ingr23_rp$roc_dat, m_mp3_rp$roc_dat, m_pb_rp$roc_dat, m_allothers_rp$roc_dat, m_noigos_rp$roc_dat)
ROC_labeldat <- rbind(m1_rp$area_under_roc, m_ingr23_rp$area_under_roc, m_mp3_rp$area_under_roc, m_pb_rp$area_under_roc, m_allothers_rp$area_under_roc, m_noigos_rp$area_under_roc)

ROC.p <- ggplot(data = ROC_dat, aes(x = x, y = y)) + geom_line() + ylim(c(0, 1)) + facet_wrap(~ Model, ncol = 1) + geom_text(data = ROC_labeldat, aes(label = paste("Area under\ncurve:",round(auc, digits = 2))), x = 1, y = 0.2, hjust = 1) + labs(title = "ROC curves", hjust = 0.5) + xlab("1 - Specificity") + ylab("Sensitivity") + theme_jk()

library(gridExtra)
pdf(file = "Output_Tables-and-Figures/crises_pr-roc.pdf", width = 6, height = 12)
grid.arrange(PRC.p, ROC.p, ncol = 2)
dev.off()

# Posterior predictive checks

m1_stan_object <- m1_roc_rstan
m1_pred_sim <- as.data.frame(m1_stan_object, pars = c("pred"))
m1_pred_sim <- reshape2::melt(m1_pred_sim)
m1_pred_sim$obs <- as.numeric(gsub("[^0-9]","",as.character(m1_pred_sim$variable)))

m1_obs_igo <- data.frame(y_obs = m1_mf[, "viol.majclash"], 
                      x_obs = m1_mf[, "igo_lev3_count_use"],
                      obs = 1:(nrow(m1_mf)))

m1_median_obs <- summarize(group_by(m1_obs_igo, x_obs),
                        median_y = mean(y_obs))

m1_pp_df <- merge(x = m1_pred_sim[, c("value", "obs")], y = m1_obs_igo, by = "obs")

ppc.p <- ggplot(data = m1_pp_df, aes(x = x_obs, y = value, group = x_obs)) + geom_boxplot(outlier.size = NA, size = 0.25) + scale_x_continuous(breaks = unique(m1_obs_igo$x_obs)) + geom_point(data = m1_median_obs, aes(x = x_obs, y = median_y), color = "red") + theme_jk() + xlab("Joint memberships in IGOs with high leverage") + ylab("Estimated probability of major clashes or war") 

ggsave(ppc.p, file = "Output_Tables-and-Figures/crises_ppc.pdf", width = 6, height = 4)

## Models with other controls

# igo_ingr23_use
m2_eq <- viol.majclash ~ igo_lev3_count_use + igo_ingr23_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally
m2_mf <- model.frame(m2_eq, data = dat)
m2 <- glm(m2_eq, data = m2_mf, family = binomial(link = "logit"))

m2_stanglm <- stan_glm(formula = m2_eq,
                       data = (m2_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m2_stanglm, file = "Output_MCMC/crises_m2_stanglm.RDS")

# Table
m2_tab <- m2_stanglm$stan_summary[, c(1, 3)]
m2_tab <- data.frame(m2_tab[c(2:9, 1, 11), ])
m2_tab$varname <- c("IGOs with high leverage",
                    "Structured IGOs", 
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m2_tab) <- c("Mean", "Standard deviation", "Variable")
m2_tab <- m2_tab[, c(3, 1, 2)]

m2_xtab <- xtable(m2_tab, digits = 2)
print(m2_xtab, file = "Output_Tables-and-Figures/crises_m2_table.tex", include.rownames = FALSE)

# with HLIGOs explicitly removed
m2b_eq <- viol.majclash ~ igo_lev3_count_use + I(igo_ingr23_use - igo_lev3_count_use) + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally
m2b_mf <- model.frame(m2b_eq, data = dat)
m2b <- glm(viol.majclash ~ ., data = m2b_mf, family = binomial(link = "logit"))

texreg(list(m2b),
       stars = 0.1,
       file = "Output_Tables-and-Figures/crises_m2b_mle_table.tex",
       custom.note = "* p < 0.05, one-tailed tests.",
       digits = 2,
       caption = "Determinants of major clashes or war in crises: comparing IGOs with high leverage to Structured IGOs. Logistic regression estimates (fit with maximum likelihood)",
       caption.above = TRUE,
       custom.coef.names = c("Intercept",
                             "IGOs with high leverage",
                             "Structured IGOs (without IGOs w/ high leverage)",
                             "Existential threat", 
                             "Territorial dispute", 
                             "Joint democracy", 
                             "Strategic rivalry",
                             "UNGA ideal point difference",
                             "Allies"))

# igo_mp3_use
m3_eq <- viol.majclash ~ igo_lev3_count_use + igo_mp3_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally
m3_mf <- model.frame(m3_eq, data = dat)
m3 <- glm(m3_eq, data = m3_mf, family = binomial(link = "logit"))

m3_stanglm <- stan_glm(formula = m3_eq,
                       data = (m3_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m3_stanglm, file = "Output_MCMC/crises_m3_stanglm.RDS")

# Table
m3_tab <- m3_stanglm$stan_summary[, c(1, 3)]
m3_tab <- data.frame(m3_tab[c(2:9, 1, 11), ])
m3_tab$varname <- c("IGOs with high leverage",
                    "Highly structured Security IGOs", 
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m3_tab) <- c("Mean", "Standard deviation", "Variable")
m3_tab <- m3_tab[, c(3, 1, 2)]

m3_xtab <- xtable(m3_tab, digits = 2)
print(m3_xtab, file = "Output_Tables-and-Figures/crises_m3_table.tex", include.rownames = FALSE)

# igo_pb_use
m4_eq <- viol.majclash ~ igo_lev3_count_use + igo_pb_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally
m4_mf <- model.frame(m4_eq, data = dat)
m4 <- glm(m4_eq, data = m4_mf, family = binomial(link = "logit"))

m4_stanglm <- stan_glm(formula = m4_eq,
                       data = (m4_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m4_stanglm, file = "Output_MCMC/crises_m4_stanglm.RDS")

# Table
m4_tab <- m4_stanglm$stan_summary[, c(1, 3)]
m4_tab <- data.frame(m4_tab[c(2:9, 1, 11), ])
m4_tab$varname <- c("IGOs with high leverage",
                    "Peace-brokering IGOs", 
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m4_tab) <- c("Mean", "Standard deviation", "Variable")
m4_tab <- m4_tab[, c(3, 1, 2)]

m4_xtab <- xtable(m4_tab, digits = 2)
print(m4_xtab, file = "Output_Tables-and-Figures/crises_m4_table.tex", include.rownames = FALSE)

# igo_allothers
m5_eq <- viol.majclash ~ igo_lev3_count_use + igo_allothers + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally
m5_mf <- model.frame(m5_eq, data = dat)
m5 <- glm(m5_eq, data = m5_mf, family = binomial(link = "logit"))

m5_stanglm <- stan_glm(formula = m5_eq,
                       data = (m5_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m5_stanglm, file = "Output_MCMC/crises_m5_stanglm.RDS")

# Table
m5_tab <- m5_stanglm$stan_summary[, c(1, 3)]
m5_tab <- data.frame(m5_tab[c(2:9, 1, 11), ])
m5_tab$varname <- c("IGOs with high leverage",
                    "All other IGOs", 
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m5_tab) <- c("Mean", "Standard deviation", "Variable")
m5_tab <- m5_tab[, c(3, 1, 2)]

m5_xtab <- xtable(m5_tab, digits = 2)
print(m5_xtab, file = "Output_Tables-and-Figures/crises_m5_table.tex", include.rownames = FALSE)

# cincdif
m6_eq <- viol.majclash ~ igo_lev3_count_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally + cincdif
m6_mf <- model.frame(m6_eq, data = dat)
m6 <- glm(m6_eq, data = m6_mf, family = binomial(link = "logit"))

m6_stanglm <- stan_glm(formula = m6_eq,
                       data = (m6_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m6_stanglm, file = "Output_MCMC/crises_m6_stanglm.RDS")

# Table
m6_tab <- m6_stanglm$stan_summary[, c(1, 3)]
m6_tab <- data.frame(m6_tab[c(2:9, 1, 11), ])
m6_tab$varname <- c("IGOs with high leverage",
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Power differential",
                    "Intercept",
                    "Log-posterior density")
names(m6_tab) <- c("Mean", "Standard deviation", "Variable")
m6_tab <- m6_tab[, c(3, 1, 2)]

m6_xtab <- xtable(m6_tab, digits = 2)
print(m6_xtab, file = "Output/crises_m6_table.tex", include.rownames = FALSE)

# tradedep_low_ln
m7_eq <- viol.majclash ~ igo_lev3_count_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally + tradedep_low_ln
m7_mf <- model.frame(m7_eq, data = dat)
m7 <- glm(m7_eq, data = m7_mf, family = binomial(link = "logit"))

m7_stanglm <- stan_glm(formula = m7_eq,
                       data = (m7_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m7_stanglm, file = "Output_MCMC/crises_m7_stanglm.RDS")

# Table
m7_tab <- m7_stanglm$stan_summary[, c(1, 3)]
m7_tab <- data.frame(m7_tab[c(2:9, 1, 11), ])
m7_tab$varname <- c("IGOs with high leverage",
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Trade dependence (lower)",
                    "Intercept",
                    "Log-posterior density")
names(m7_tab) <- c("Mean", "Standard deviation", "Variable")
m7_tab <- m7_tab[, c(3, 1, 2)]

m7_xtab <- xtable(m7_tab, digits = 2)
print(m7_xtab, file = "Output_Tables-and-Figures/crises_m7_table.tex", include.rownames = FALSE)

# gdppc_low_ln
m8_eq <- viol.majclash ~ igo_lev3_count_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally + gdppc_low_ln
m8_mf <- model.frame(m8_eq, data = dat)
m8 <- glm(m8_eq, data = m8_mf, family = binomial(link = "logit"))

m8_stanglm <- stan_glm(formula = m8_eq,
                       data = (m8_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m8_stanglm, file = "Output_MCMC/crises_m8_stanglm.RDS")

# Table
m8_tab <- m8_stanglm$stan_summary[, c(1, 3)]
m8_tab <- data.frame(m8_tab[c(2:9, 1, 11), ])
m8_tab$varname <- c("IGOs with high leverage",
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "GDP per capita (lower)",
                    "Intercept",
                    "Log-posterior density")
names(m8_tab) <- c("Mean", "Standard deviation", "Variable")
m8_tab <- m8_tab[, c(3, 1, 2)]

m8_xtab <- xtable(m8_tab, digits = 2)
print(m8_xtab, file = "Output_Tables-and-Figures/crises_m8_table.tex", include.rownames = FALSE)

# continentFactor
m9_eq <- viol.majclash ~ igo_lev3_count_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally + continentFactor
m9_mf <- model.frame(m9_eq, data = dat)
m9 <- glm(m9_eq, data = m9_mf, family = binomial(link = "logit"))

m9_stanglm <- stan_glm(formula = m9_eq,
                       data = (m9_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m9_stanglm, file = "Output_MCMC/crises_m9_stanglm.RDS")

# Table
m9_tab <- m9_stanglm$stan_summary[, c(1, 3)]
m9_tab <- data.frame(m9_tab[c(2:12, 1, 13), ])
m9_tab$varname <- c("IGOs with high leverage",
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Africa (vs. different continents)",
                    "Americas (vs. different continents)",
                    "Asia/Oceania (vs. different continents)",
                    "Europe (vs. different continents)",
                    "Intercept",
                    "Log-posterior density")
names(m9_tab) <- c("Mean", "Standard deviation", "Variable")
m9_tab <- m9_tab[, c(3, 1, 2)]

m9_xtab <- xtable(m9_tab, digits = 2)
print(m9_xtab, file = "Output_Tables-and-Figures/crises_m9_table.tex", include.rownames = FALSE)

# coldWar
m10_eq <- viol.majclash ~ igo_lev3_count_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally + coldWar
m10_mf <- model.frame(m10_eq, data = dat)
m10 <- glm(m10_eq, data = m10_mf, family = binomial(link = "logit"))

m10_stanglm <- stan_glm(formula = m10_eq,
                        data = (m10_mf),
                        family = binomial(link = "logit"),
                        prior = normal(location = 0, scale = 10), 
                        prior_intercept = normal(0, 10),
                        seed = 20170207,
                        cores = 4)

saveRDS(m10_stanglm, file = "Output_MCMC/crises_m10_stanglm.RDS")

# Table
m10_tab <- m10_stanglm$stan_summary[, c(1, 3)]
m10_tab <- data.frame(m10_tab[c(2:9, 1, 11), ])
m10_tab$varname <- c("IGOs with high leverage",
                     "Existential threat", 
                     "Territorial dispute", 
                     "Joint democracy", 
                     "Strategic rivalry",
                     "UNGA ideal point difference",
                     "Allies",
                     "Cold War",
                     "Intercept",
                     "Log-posterior density")
names(m10_tab) <- c("Mean", "Standard deviation", "Variable")
m10_tab <- m10_tab[, c(3, 1, 2)]

m10_xtab <- xtable(m10_tab, digits = 2)
print(m10_xtab, file = "Output_Tables-and-Figures/crises_m10_table.tex", include.rownames = FALSE)

# Random effects for time
m14_eq <- viol.majclash ~ (igo_lev3_count_use | decade) + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally
m14_mf <- model.frame(viol.majclash ~ igo_lev3_count_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally + decade, data = dat)
m14 <- glmer(m14_eq, data = dat, family = binomial(link = "logit"))

m14_stanglm <- stan_glmer(formula = m14_eq,
                          data = (dat),
                          family = binomial(link = "logit"),
                          prior = normal(location = 0, scale = 10), 
                          prior_intercept = normal(0, 10),
                          seed = 20170207,
                          cores = 4)

saveRDS(m14_stanglm, file = "Output_MCMC/crises_m14_stanglm.RDS")

m14_mcmc <- as.data.frame(m14_stanglm$stanfit)
m14_re_mcmc <- m14_mcmc[ , grepl( "igo_lev3_count_use decade:" , names(m14_mcmc) ) ]
m14_re_mcmc$`b[igo_lev3_count_use decade:_NEW_decade]` <- NULL
table(m14_mf$decade)
# 0   1   2   3   4   5   6 
# 17  51  81  98 101 136  42 
names(m14_re_mcmc) <- c("1947-1949 (17 crises)", "1950-1959 (51 crises)", "1960-1969 (81 crises)", "1970-1979 (98 crises)", "1980-1989 (101 crises)", "1990-2005 (178 crises)")

m14_re_mcmc_probs <- data.frame(variable = names(m14_re_mcmc), 
                                perc_below_0 = apply(m14_re_mcmc, MARGIN = 2, FUN = function(x) sum(ifelse(x < 0, 1, 0)) / length(x)),
                                perc_larger = apply(m14_re_mcmc, MARGIN = 2, FUN = function(x) ifelse(sum(ifelse(x < 0, 1, 0)) / length(x) > 0.5, sum(ifelse(x < 0, 1, 0)) / length(x), sum(ifelse(x > 0, 1, 0)) / length(x))),
                                perc_below_0_label = paste(round(apply(m14_re_mcmc, MARGIN = 2, FUN = function(x) sum(ifelse(x < 0, 1, 0)) / length(x)) * 100, digits = 1), "% < 0", sep = ""), 
                                xpos = apply(m14_re_mcmc, MARGIN = 2, FUN = function(x) quantile(x, probs = 0.01)))

m14_re_mcmc_long <- reshape2::melt(m14_re_mcmc)
p <- ggplot(data = m14_re_mcmc_long, aes(x = value)) + geom_density(color = NA, fill = "lightgray") + geom_text(data = m14_re_mcmc_probs, aes(x = xpos, y = 0.1, label = perc_below_0_label), vjust = 0) + facet_wrap(~ variable, scales = "free", ncol = 1) + geom_vline(xintercept = 0) + theme_jk() + xlab("Estimate for decade-specific logit coefficients on IGOs with high leverage") + ylab("") + scale_y_continuous(breaks = NULL)

ggsave(p, file = "Output_Tables-and-Figures/crises_hligo_re.pdf", width = 6, height = 8)

# Model comparison: are the REs worth it? 
loo_m14 <- loo(m14_stanglm, cores = 4)
compare(loo_m1, loo_m14)
# elpd_diff  se 
# 36.6       9.0

# Interact HLIGOs and commitment problem (power differential)
m15_eq <- viol.majclash ~ igo_lev3_count_use * cincdif + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally
m15_mf <- model.frame(m15_eq, data = dat)
m15 <- glm(m15_eq, data = m15_mf, family = binomial(link = "logit"))

m15_mf$cincdif_q3 <- cut(m15_mf$cincdif, breaks = quantile(m15_mf$cincdif, probs = c(0, 0.33, 0.66, 1)), include.lowest = TRUE)
levels(m15_mf$cincdif_q3) <- c("[0, 0.002]", "(0.002, 0.01]", "(0.01, 0.3]")

m15b_eq <- viol.majclash ~ igo_lev3_count_use * cincdif_q3 + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally
m15b <- glm(m15b_eq, data = m15_mf, family = binomial(link = "logit"))

m15_stanglm <- stan_glm(formula = m15b_eq,
                        data = (m15_mf),
                        family = binomial(link = "logit"),
                        prior = normal(location = 0, scale = 10), 
                        prior_intercept = normal(0, 10),
                        seed = 20170207,
                        cores = 4)

saveRDS(m15_stanglm, file = "Output_MCMC/crises_m15_stanglm.RDS")

m15_mcmc <- as.data.frame(m15_stanglm$stanfit)
igo_lowcincdif <- data.frame(x = m15_mcmc[, "igo_lev3_count_use"], var = "No or small CINC difference (0 to 0.002)")
igo_midcincdif <- data.frame(x = m15_mcmc[, "igo_lev3_count_use"] + m15_mcmc[, "igo_lev3_count_use:cincdif_q3(0.002, 0.01]"], var = "Medium CINC difference (0.002 to 0.01)")
igo_hicincdif <- data.frame(x = m15_mcmc[, "igo_lev3_count_use"] + m15_mcmc[, "igo_lev3_count_use:cincdif_q3(0.01, 0.3]"], var = "High CINC difference (0.01 to 0.3)")
m15_interact_mcmc <- rbind(igo_lowcincdif, igo_midcincdif, igo_hicincdif)

perc_below_0 <- c(sum(ifelse(igo_lowcincdif$x < 0, 1, 0)) / length(igo_lowcincdif$x), sum(ifelse(igo_midcincdif$x < 0, 1, 0)) / length(igo_midcincdif$x), sum(ifelse(igo_hicincdif$x < 0, 1, 0)) / length(igo_hicincdif$x))

m15_interact_mcmc_probs <- data.frame(var = c("No or small CINC difference (0 to 0.002)", "Medium CINC difference (0.002 to 0.01)", "High CINC difference (0.01 to 0.3)"), 
                                      perc_below_0 = c(sum(ifelse(igo_lowcincdif$x < 0, 1, 0)) / length(igo_lowcincdif$x), sum(ifelse(igo_midcincdif$x < 0, 1, 0)) / length(igo_midcincdif$x), sum(ifelse(igo_hicincdif$x < 0, 1, 0)) / length(igo_hicincdif$x)),
                                      perc_below_0_label = paste(round(perc_below_0 * 100, digits = 1), "% < 0", sep = ""),
                                      xpos = c(-0.8, -0.6, -0.9))

p <- ggplot(data = m15_interact_mcmc, aes(x = x)) + geom_density(color = NA, fill = "lightgray") + geom_text(data = m15_interact_mcmc_probs, aes(x = xpos, y = 0.6, label = perc_below_0_label), hjust = 0) + facet_wrap(~ var, scales = "free", ncol = 1) + geom_vline(xintercept = 0) + theme_jk() + xlab("Estimate for conditional effect (dy/dx) of IGOs with high leverage") + ylab("") + scale_y_continuous(breaks = NULL)

ggsave(p, file = "Output_Tables-and-Figures/crises_igo-cinc_interact.pdf", width = 5, height = 8)

## Table of MLE results

texreg(list(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10),
        stars = 0.1,
        file = "Output_Tables-and-Figures/crises_mle-results.tex",
        # reorder.coef = c(),
        custom.note = "* p < 0.05, one-tailed tests.",
        digits = 2,
        caption = "Determinants of major clashes or wars during crises: logistic regression estimates (fit with maximum likelihood)",
        caption.above = TRUE,
        custom.coef.names = c("Intercept",
                              "IGOs with high leverage",
                              "Existential threat", 
                              "Territorial dispute", 
                              "Joint democracy", 
                              "Strategic rivalry",
                              "UNGA ideal point difference",
                              "Allies",
                              "Structured IGOs",
                              "Highly structured Security IGOs", 
                              "Peace-brokering IGOs",
                              "All other IGOs",
                              "Power differential",
                              "Trade dependence (lower)",
                              "GDP per capita (lower)",
                              "Africa (vs. different continents)",
                              "Americas (vs. different continents)",
                              "Asia/Oceania (vs. different continents)",
                              "Europe (vs. different continents)",
                              "Cold war"))

#############################################
# 2. Bayesian Model Averaging
#############################################

bma_dat <- model.frame(viol.majclash ~ igo_lev3_count_use + igo_ingr23_use + igo_mp3_use + igo_pb_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally + cincdif, data = dat)
names(bma_dat) <- c("Major clashes or war",
                    "IGOs with high leverage",
                    "Structured IGOs",
                    "Highly structured Security IGOs", 
                    "Peace-brokering IGOs",
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Power differential")

crises_bma <- bic.glm(y = bma_dat[, 1], x = bma_dat[, c(2:12)], 
                      glm.family = binomial(link = "logit"))

pdf("Output_Tables-and-Figures/crises_bma.pdf", width = 12, height = 8)
plotJK.bic.glm(crises_bma, mfrow = c(3, 4))
dev.off()

#############################################
# 3. IGO characteristics
#############################################

char_dat <- model.frame(viol.majclash ~ igo_lev3_count_use + igo_ingr23_use + igo_mp3_use + igo_pb_use + highgrav + territory + jointpol7 + jointpol6 + rivalry_th + absidealdiff + atopally + cincdif, data = dat)

char_dat$idealpointdiff_std <- scales::rescale(x = char_dat$absidealdiff, to = c(0, 1))
char_dat$cincdif_std <- scales::rescale(x = char_dat$cincdif, to = c(0, 1))

char_dat$hligo_quartiles <- cut(char_dat$igo_lev3_count_use, breaks = c(0, 1, 2, 3, 7), include.lowest = TRUE)
char_dat$hligo_quartiles <- factor(char_dat$hligo_quartiles, labels = c("1st", "2nd", "3rd", "4th"))

char_dat$ingr23_quartiles <- cut(char_dat$igo_ingr23_use, breaks = quantile(char_dat$igo_ingr23_use, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE), include.lowest = TRUE)
char_dat$ingr23_quartiles <- factor(char_dat$ingr23_quartiles, labels = c("1st", "2nd", "3rd", "4th"))

char_dat$mp3_quartiles <- factor(char_dat$igo_mp3_use)
char_dat$mp3_quartiles <- factor(char_dat$mp3_quartiles, labels = c("1st", "2nd", "3rd", "4th"))

char_dat$pb_quartiles <- ifelse(char_dat$igo_pb_use <= 1, 1, 
                                ifelse(char_dat$igo_pb_use == 2, 2,
                                       ifelse(char_dat$igo_pb_use == 3, 3,
                                              ifelse(char_dat$igo_pb_use >= 4, 4, NA))))
char_dat$pb_quartiles <- factor(char_dat$pb_quartiles, labels = c("1st", "2nd", "3rd", "4th"))

hligo_sum <- summarize(group_by(char_dat[is.na(char_dat$hligo_quartiles) == FALSE, ], hligo_quartiles), 
                       obs_ratio = length(atopally) / nrow(char_dat[is.na(char_dat$hligo_quartiles) == FALSE, ]), 
                       obs_count = length(atopally),
                       democracies_mean = mean(jointpol6, na.rm = TRUE), 
                       rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                       idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                       allies_mean = mean(atopally, na.rm = TRUE), 
                       cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(hligo_sum) <- c("Joint IGO memberships", "% of observations", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
hligo_sum$`IGO type` <- 1 # "IGOs with high leverage"

# igo_ingr23_use
ingr23_sum <- summarize(group_by(char_dat[is.na(char_dat$ingr23_quartiles) == FALSE, ], ingr23_quartiles), 
                        obs_ratio = length(atopally) / nrow(char_dat[is.na(char_dat$ingr23_quartiles) == FALSE, ]), 
                        obs_count = length(atopally),
                        democracies_mean = mean(jointpol6, na.rm = TRUE), 
                        rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                        idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                        allies_mean = mean(atopally, na.rm = TRUE), 
                        cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(ingr23_sum) <- c("Joint IGO memberships", "% of observations", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
ingr23_sum$`IGO type` <- 2 # "Structured IGOs"

# igo_mp3_use
mp3_sum <- summarize(group_by(char_dat[is.na(char_dat$mp3_quartiles) == FALSE, ], mp3_quartiles), 
                     obs_ratio = length(atopally) / nrow(char_dat[is.na(char_dat$mp3_quartiles) == FALSE, ]), 
                     obs_count = length(atopally),
                     democracies_mean = mean(jointpol6, na.rm = TRUE), 
                     rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                     idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                     allies_mean = mean(atopally, na.rm = TRUE), 
                     cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(mp3_sum) <- c("Joint IGO memberships", "% of observations", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
mp3_sum$`IGO type` <- 3 # "Highly structured Security IGOs"

# igo_pb_use
pb_sum <- summarize(group_by(char_dat[is.na(char_dat$pb_quartiles) == FALSE, ], pb_quartiles), 
                    obs_ratio = length(atopally) / nrow(char_dat[is.na(char_dat$pb_quartiles) == FALSE, ]), 
                    obs_count = length(atopally),
                    democracies_mean = mean(jointpol6, na.rm = TRUE), 
                    rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                    idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                    allies_mean = mean(atopally, na.rm = TRUE), 
                    cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(pb_sum) <- c("Joint IGO memberships", "% of observations", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
pb_sum$`IGO type` <- 4 # "Peace-brokering IGOs"

igo_sum <- rbind(hligo_sum, ingr23_sum, mp3_sum, pb_sum)
names(igo_sum) <- c("Joint IGO memberships", "% of observations", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential", "IGO type")
igo_sum$`IGO type` <- factor(igo_sum$`IGO type`, labels = c("IGOs with high leverage", "Structured IGOs", "Highly structured Security IGOs", "Peace-brokering IGOs"))
igo_sum <- igo_sum[is.na(igo_sum$`Joint IGO memberships`) == FALSE, ]

igo_sum$`% of observations` <- NULL
igo_sum_long <- reshape2::melt(igo_sum, id.vars = c("Joint IGO memberships", "IGO type", "Observations"))

# igo_sum$`Observations` <- NULL
# igo_sum_long <- reshape2::melt(igo_sum, id.vars = c("Joint IGO memberships", "IGO type", "% of observations"))

p <- ggplot(data = igo_sum_long, aes(x = `Joint IGO memberships`, y = value)) + facet_grid(variable ~ `IGO type`) + geom_line(aes(group = 1)) + geom_point(aes(size = Observations), shape = 21, fill = "white") + theme_jk() + xlab("Count of joint IGO memberships (quartile-based bins)") + ylab("Average (proportions or means standardized to 0-1 range)") + labs(size = "Crises\nin bin")

ggsave(p, file = "Output_Tables-and-Figures/crises_igo-profiles.pdf", width = 9.75, height = 10)

# For analyses with crises collapsed, see crises_coll_analysis.R

# End of file.